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ABSTRACT 


Econometricians  must  choose  between  many  methods  for  estimating  p,  the 
autocorrelation  coefficient,  in  a  first  order  autoregressive  process.  This  thesis  examines 
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The  autocorrelation  coeficient,  p,  was  varied  from  .2  to  .9  and  each  method  provided 
estimates  of  p  and  P,  the  regression  coefficient,  for  1000  replications.  The  results 
presented  here  are  similar  to  those  found  in  previous  comparisons.  Specifically, 
Ordinary  Least  Squares  was  found  to  be  an  efficient  estimator  of  p  when 
autocorrelation  is  present  only  to  a  slight  degree.  Of  the  four  estimators  examined,  the 
performance  of  Theil-Nagar  proved  superior  in  estimating  both  p  and  p  far  small 
values  of  the  correlation  coeficient.  Beach-MacKinnon,  on  the  other  hand,  while 
containing  a  large  bias  in  the  estimation  of  p,  is  the  more  efficient  estimator  of  p  for 
large  values  of  p. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Autocorrelation  exists  in  a  regression  model  when  the  error  terms  are  no  longer 
independent  but  are  correlated.  In  the  examination  of  time  scries  data  autocorrelation 
is  a  common  phenomenon  and  can  lead  to  problems  if  Ordinary  Least  Squares  (OLS) 
estimation  procedures  are  used.  The  purpose  of  this  thesis  is  to  examine  and  compare 
four  different  estimates  of  the  autocorrelation  coefficient,  p,  the  estimation  of  which  is 
essential  to  the  resolution  of  OLS  deficiencies.  The  four  estimators  to  be  examined  arc 
the  Durbin-Watson,  Theil-Nagar,  Beach-MacKinnon,  and  Prais-Winsten. 

B.  PROBLEM  STATEMENT 

In  the  standard  regression  model  y  =  XP  +e,  y  is  a  Txl  vector  of  observations  of 
a  dependent  variable,  X  is  a  TxK  design  matrix  and  P  is  a  Kxl  vector.  The  variable  e 
is  a  Txl  vector  of  unobservable  random  errors  with  E(e)  =  0  and  covariance  matrix, 
E(ee')=  g^Ij.  Thus,  in  the  standard  model,  the  random  vector  e  contains  elements 
which  are  pairwise  uncorrclatcd  with  identical  means  and  variances.  In  the  presence  of 
autocorrelation  this  strong  assumption  is  violated.  That  is,  the  error  terms  arc  no 
longer  independent  but  arc  correlated.  The  regression  model  becomes, 

yt  =  Xtp  +  et  t- 1,2 . T  (cqn  1.1) 

where  et  -  p  ct_j  +  vt, 

E(vt)  =  0,  and 
E(vv')  =  <t2I  . 

This  is  known  as  a  first  order  autoregressive  or  AR(1)  process.  As  illustrated  by 
equation  1.1,  ct  is  expressed  linearly  in  terms  of  the  ct_j  and  another  random  error 
term  vt.  The  assumption  of  zero  mean  and  constant  variance  provides  vt  with  all  the 
nice  properties  of  et  in  the  standard  model.  This  process  may  occur  for  a  variety  of 
reasons,  some  of  which  arc: 

1  Omitted  explanatorv  variables.  If  a  correlated  explanatory  variable  has  been 
excluded  from  the  "design  matrix  its  exclusion  will  be  reflected  in  the  correlation 
of  the  random  variable"!:. 

2  Mispec.ificalion  of  the  mathematical  form  of  the  model.  If  the  wrong  mathematical 
relationship  is  chosen  the  values  of  c  may" be  dependent. 


3  Interpolations  in  the  statistical  observations.  If  the  observational  data  is  smoothed 
autocorrelation  may  result. 

4  Mispecification  of  the  true  random  error.  Dependence  among  the  error  terms  may 
occur  naturally.  [Ref.  l:p.  204] 

Utilizing  OLS  to  estimate  the  regression  coefficient,  P,  in  the  presence  of  an  AR(1) 
process  can  lead  to  problems.  Generally,  there  are  two  consequences  to  consider.  The 
first  is  that  the  OLS  estimator  of  the  coefficients  will  be  unbiased  but  will  not  be  very 
efficient.  The  second  consequence  is  that  the  OLS  variance  estimator  is  biased.  For 
these  reasons  it  is  useful  to  investigate  other  methods  to  estimate  p  [Ref.  2:p.  439]. 

C.  ESTIMATORS 

When  p  is  known,  the  process  is  easily  accounted  for  using  Generalized  Least 
Squares  or  Weighted  Least  Squares  methods  [Ref.  3].  However,  the  usual  situation  is 
that  p  is  unknown  and  must  be  estimated.  A  number  of  methods  have  been  proposed 
to  estimate  p  and  properly  account  for  OLS  deficiencies  in  estimating  p.  Chapter  2  will 
develop  the  four  estimators  mentioned  above  and  examine  the  autocorrelation  process. 

D.  SIMULATION 

Each  of  the  estimators  considered  here  have  the  same  asymptotic  properties 
therefore  any  decision  on  which  one  to  use  must  be  based  on  small  sample  analysis  and 
Monte-Carlo  evidence.  Therefore,  a  simulation  will  be  created  in  which  the  data  is 
generated  according  to  guidelines  presented  in  previous  studies  with  equation  1.1  as  the 
model.  The  actual  values  of  p  will  be  varied  from  .2  to  .9.  The  four  estimation 
techniques  will  then  provide  estimates  of  p  and  P  for  1000  replications. 

E.  MEASURE  OF  EFFECTIVENESS 

To  provide  an  indication  of  which  estimator  performs  best  the  mean  square  error 

A  /S 

of  both  p  and  P  will  be  estimated  for  each  estimator.  Prior  results  for  different  sets  of 
estimators  indicate  that  no  one  estimator  will  prove  superior  over  the  entire  range  of  p 
but  that  one  or  two  may  out  perform  the  others  over  specific  intervals. 


II.  ESTIMATION 


A.  GENERAL 

This  chapter  attempts  to  explore  the  theory  behind  both  the  first  order  process 
and  four  estimators  developed  to  properly  account  for  it.  Three  of  these 
(Durbin- Watson,  Thicl-N'agar,  and  Prais-Winsten)  are  categorized  as  estimated 
generalized  least  squares  estimators.  The  fourth  (Beach-MacKinnon)  is  a  maximum 
likelihood  estimator. 

B.  PROCESS 


The  first  order  process  can  be  written  as 

>'t  =  XtP  +  et  t=  1,2, ...,T  (cqn  2.1) 

where  e(  =  p  c,  <  +  vt, 

E(vt)  =  0, 

E(vv')  =  <y2l, 

E(vt2)  =  <rv2,  and 
E(vtvs)~0  for  s*t  . 

The  parameter  p  is  generally  unknown  and  along  with  P  must  be  estimated.  The 
statistical  properties  of  the  random  error,  v,  listed  in  equation  2.1  arc  identical  to  those 
listed  for  c  in  the  general  linear  model.  The  statistical  properties  of  c  under  these  new 
assumptions  arc  quite  different.  Judge  [Ref.  4:p.  438]  shows  that 


E(et)  =  £  piE(vt_i)  =  0 


(eqn  2.2) 


E(et2)=<Te2  =  <Tv2/(l-p)2  . 


(cqn  2.3) 


The  covariance  between  errors  s  periods  apart  is  no  longer  zero  and  is  given  by 


E(etct_s)=  E(c,  +  c.)  =  (ps<r2  )'( 1-p2)  . 


(cqn  2.4) 


The  covariance  matrix  for  e  is  now  easily  written 

as 

O  =  E(ec')  = 

1 

p  .  .  . 

pt-‘ 

(eqn  2.5) 

P 

1  . 

pT-2 

=  cr2v/(!-p2) 

T-l 

P 

T-2 

P  •  •  • 

1 

or  utilizing  the  iollowing  convention, 

O  =  <tv2,I' 

(eqn  2.6) 

where  ’I'  = 

1 

a 

P  P^ 

•  P1-1 

P 

1  p2  . 

•  PT'2 

P2 

P  1 

•  PT-3 

=  l/d -p2) 

pT-1 

p'l  -2  pT-3 

1 

• 

Thus,  the  assumptions  made  about  the  error  term,  e,  in  the  standard  linear  model 
no  longer  hold  for  the  autoregressive  case.  Specifically,  due  to  autocorrelation  the  error 
covariance  matrix  is  no  longer  written  as  a"!  but  is  now  <y  yT  . 

When  an  attempt  is  made  to  perform  a  least  squares  fit  to  the  data  in  the  presence 

of  an  AR(!)  process  there  arc  two  problems  to  consider. 

A  1 

1  The  least  squares  estimator  H  =  (X'Xf'X'v  will  be  unbiased  but  will  not  be  very 
efficient. 

■A")  I  f 

2  The  least  squares  covariance  matrix  <r~(XA«j  with  (T‘  =  (y-Xb)'(y-Xb)  (T-k) 
will  be  a  biased  estimator  ul  the  variance  ol  p. 


In  the  presence  of  positive  autocorrelation  Judge  |Ref.  4:p.  439]  notes  that  with 


OLS  estimation  the  bias  of  the  standard  error  of  p  will  very  likely  appear  as  an 


II 


underestimate.  Park  and  Mitchell  [Ref.  5.p.  16]  warn  that  OLS  seriously 
underestimates  the  variance  of  P  for  p  >  0.4.  This  understatement  makes  the 
estimates  themselves  appear  much  more  significant  than  they  actually  are  and  makes 
hypothesis  testing  of  the  slope  coefficients  unreliable. 

C.  METHODS  OF  ESTIMATION 

1.  Generalized  Least  Squares  Estimation 

When  apriori  information  is  available  about  *F,  the  most  convenient  estimate 

A. 

for  the  regression  coefficient,  p,  is  obtained  by  applying  least  squares  estimation 
techniques  to  the  transformed  model, 

Y*  =  X*p  +  e*  (eqn  2.7) 

where  Y*  =  PY 
X*  =  PX 
e*  =  Pe. 

The  transformation  matrix  P  is  the  TxT  matrix 


V  1-p2  000.  .  .  0 

-p  1  0  0  .  .  0 

0  -p  1  0  .  .  0 

P  =  0  0  -p  1  ...  0 

000.  .  .  -pl 

where  P'P  = 


This  method  is  known  as  the  Generalized  Least  Squares  (GLS)  estimation. 

2.  Estimated  Generalized  Least  Squares 

The  usual  case  is  that  p  is  unknown  and  must  be  estimated.  Once  an  estimate 
-a  A 

for  p  (p)  is  computed  one  can  substitute  p  into  the  P  matrix  and  proceed  with  the  GLS 
method  outlined  above.  This  is  known  as  Estimated  Generalized  Least  Squares 
(EGLS)  estimation.  The  computational  form  of  the  alternative  estimators  for  p 
discussed  are  as  follows: 


a.  Durbin-  Watson 
The  statistic 


d  ■  Z  (vct.i)2/S V  t=1 . T 

wliere  et  = 


(eqn  2.8) 


is  often  used  to  test  for  first  order  autoregressive  errors.  As  the  number  of 
observations  (T)  increases  it  can  be  demonstrated  that  d  approaches  the  least  squares 
estimator  of  p  or 

^  =  1-  (d/2)  .  (cqn  2.9) 


The  Durbin-Watson  statistic  is  provided  by  most  least  squares  computer  packages  and 
is  very  easy  to  use.  It  also  is  an  example  of  a  two-stage  estimator.  That  is,  it  first 
estimates  the  correlation  parameter  and  then  uses  this  estimate  to  compute  the 
generalized  least  squares  estimates  for  p. 

b.  Theil-Nagar 

A  modification  of  the  Durbin-Watson  estimator  suggested  by  Henri 
Theil  and  A.  L.  Nagar  is 


P  =  (T2(l-(d/2))  +  K2)  /(T2  -  K2).  (cqn  2.10) 

Theil  and  Nagar  claim  that  this  estimator  is  an  improvement  over  Durbin-Watson  if 
the  first  and  second  differences  of  the  explanatory  variables  arc  small  when  compared 
to  their  corresponding  ranges  [Ref.  6).  Like  Durbin-Watson,  it  also  is  a  two-stage 
estimator. 

c.  Prais-Winslen 

A  minimum  sum  of  squares  approach  to  estimating  p  yields, 


A 


7- , 


P  =  I  Vt-l/I  e't2  t=l . T 


I'J  .  t  =  ‘  A, 

where  c~t  =  yt-XtP  . 


(cqn  2.1 1) 
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This  estimator  can  be  employed  in  both  a  two  step  and  an  iterative  procedure.  This 
paper,  however,  considers  only  the  following  iterative  form: 

1.  Set"?  =  0. 

2.  Transform  the  variables  in  accordance  with  the  transformation  matrix  and 
equation  2.7. 

3.  Calculate  the  least  squares  estimate  of  P  conditional  on  p. 

4.  Calculate  the  estimate  of  p  conditional  on  P  by  using  equation  2.1 1. 

5.  If  the  absolute  difference  in'?"  from  the  previous  iteration  is  sufficiently  small  (less 
than  0.00001)  stop.  If  not  go  to  step  2.  [Ref.  7:p.  2J 


3.  Maximum  Likelihood  Estimation 

A  maximum  likelihood  (ML)  estimator  is  the  value  of  0  which  maximizes  the 
value  of  the  likelihood  function  L(G).  Under  the  assumption  that  Y  has  a  multivariate 
normal  distribution  with  mean  XP  and  covariance  matrix  <r2,P,  the  likelihood  function 
is 


L(p,p,<r2)=  C-  (l/2<r2v)(y-XP)'T-1(y-XP)  (cqn  2.12) 

where  C  =  -(T/2)ln<r2y  +  (l/2)In  (1-p2) . 

The  ML  estimators  for  P,  p,  and  <ry2  are  those  values  for  which, 

5L/dp  =  0  ,  dL/dp  =  0  ,  H/5  u2y  =  0  .  (eqn  2.13) 


Solutions  to  equations  2.13  are  very  difficult  to  derive.  Beach  and  MacKinnon 
[Ref.  8:p.  54|  use  an  ML  estimator  for  <5y2  and  substitute  into  equation  2.12.  The 
result  is  the  concentrated  likelihood  function, 

L(P,p)=  K-(T/2)ln((y-XP)'T',(y-XP)(l-p2)1/T)  (cqn  2.14) 

where  K  =  (T/2)ln(T)-(T/2) . 

They  suggest  maximizing  L(P,p)  with  respect  to  P  with  p  held  constant  and  then  to 
maximize  with  respect  to  p  with  P  held  constant.  An  algorithm  to  derive  this  ML 


estimate  is 


2.  Transform  the  variables  in  accordance  with  equation  2.7. 

3.  Calculate  the  least  squares  estimate  of  P  conditional  on  pi 

4.  Calculate  the  ML  estimate  of  p  conditional  on  P  by  solving  a  cubic  equation  of 
the  untransformed  residuals,  (see  [Ref.  8j  for  details) 

5.  If  the  absolute  difference  in'p'  from  the  previous  iteration  is  sufficiently  small  (less 
than  0.00001)  stop.  If  not,  go  to  step  2  [Ref.  7).  (Note:  The  same  procedure  was 
employed  for  iterative  Prais-Winsten  method  except  that  equation  2. 1 1  was  used  to 
estimate  p.) 

This  is  not  a  comprehensive  listing  of  all  available  estimators  for  a  first  order 
process.  Other  estimators  are  listed  in  Judge  [Ref.  41. 


III.  COMPARISON 


A.  GENERAL 

The  finite  sampling  properties  of  the  estimators  listed  here  have  not  been  derived. 
Choice  of  which  estimator  to  use  might  be  based  on  evidence  obtained  from  Monte 
Carlo  simulations.  This  chapter  explains  a  simulation  used  and  provides  a  synopsis  of 
comparisons  reported  in  the  literature. 

B.  PREVIOUS  COMPARISONS 

There  have  been  a  number  of  studies  of  estimators  for  p  .  Each  has  concluded 
that  OLS  has  serious  deficiencies  in  the  presence  of  autocorrelation.  The  majority  of 
these  papers  have  settled  on  two  points.  First,  particularly  in  small  sample  sizes 
(T  <  50)  it  is  best  to  use  estimators  that  consider  all  T  observations.  Rao  and 
Grilitches  concluded  that  using  estimators  such  as  Cochrane-Orcutt  that  ignore  the 
first  observation  can  lead  to  a  substantial  loss  of  efficiency  [Ref.  9:p.  269],  These 
results  were  further  substantiated  by  Beach  and  MacKinnon.  In  an  attempt  to  develop 
a  computationally  efficient  algorithm  to  maximize  the  likelihood  function  they 
discovered  (for  p=  0.6,  0.8,  0.99)  significant  gains  in  efficiency  to  be  made  by 
employing  the  first  observation.  Some  of  these  gains  arc  in  the  neighborhood  of  700 
percent  [Ref.  8:p.  55).  Park  and  Mitchell  concluded  that  retention  of  this  first 
observation  substantially  reduces  the  risk  of  collinearity  as  p  approaches  0.9  [Ref  5:p. 
10[.  Kobayashi  verified  theoretically  the  experimental  results  of  Park  and  Mitchell.  By 
computing  the  asymptotic  variances  of  several  estimators  he  demonstrated  that  the  loss 
of  efficiency  of  the  Cochrane-Orcutt  method  was  due  primarily  to  ignoring  the  first 
observation.  [Ref.  10:p.  951]. 

The  second  point  is  that  the  Prais-Winstcn  solution  techniques  outperform  many- 
comparable  estimators  of  the  correlation  parameter.  Spitzer  concluded  that 
Prais-Winstcn  "appeared  to  be  the  best  of  all  the  two  stage  estimators."  [Ref.  ll:p. 
44],  Park  and  Mitchell  in  a  later  study  comparing  Beach- MacKinnon  with  the  iterative 
Prais-Winstcn  estimator  concluded  that  the  iterative  Prais-Winstcn  performs 
"appreciably  better  in  estimating  the  autocorrelation  coefficient  p"  [Ref  7:p.  5]. 
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Although  there  were  no  studies  found  specifically  comparing  the  four  estimators 
presented  here,  each  has  demonstrated  a  superiority  to  OLS  in  the  presence  of  a  first 
order  process. 

C.  MODEL  AND  DATA  GENERATION 

Equation  2.1  was  utilized  as  the  model  with  the  first  term  in  the  vector  e 
generated  in  the  following  fashion, 

el=  Vj/Cl-p)1''2.  (eqn  3.1) 

In  order  to  conform  with  previous  comparisons,  the  data  utilized  in  this 
experiment  is  identical  to  that  used  in  Beach  and  MacKinnon  [Ref.  SJ.  Two  sample 
sizes  of  20  and  50  observations  were  used.  The  untrended  explanatory  variable,  X,  was 
drawn  from  N(0,  0.0625)  and  the  random  error,  vt,  was  drawn  from  N(0,  0.0036). 
Although  autocorrelation  in  theory  may  be  positive  or  negative,  in  econometric  data  it 
is  almost  always  positive  [Ref.  l:p.  201],  For  this  reason  p  was  varied  from  0.2  to  0.9. 

D.  VALIDATION 

The  data  generation  program  was  checked  to  ensure  the  normality  of  c  using  the 
Chi  Square  Goodness  of  Fit  test.  The  normality  assumption  was  accepted  at  a  0.36S4 
level.  Finally,  in  order  to  ensure  each  estimator  performed  properly  the  random 
portion  of  the  model,  specifically  the  random  variable  V,  was  removed.  This  allowed 
the  estimators  to  function  in  a  deterministic  fashion.  Data  were  then  generated  and 
submitted  to  each  estimator  for  values  of  p  equal  0.2,  0.6,  0.8.  The  results  are 
presented  in  Table  I  ,  illustrate  that  the  estimators  are  functioning  properly. 

E.  SIMULATION 

For  each  run  the  values  of  the  regression  coefficients,  Pq  and  Pj,  were  set  to  1 
and  1.  The  variables  Xt  and  vt  were  drawn  from  the  normal  distributions  discussed 
earlier.  The  dependent  variable  yt  was  calculated  using  equation  2.1.  Since  the 
ultimate  objective  was  to  generate  residuals  to  send  to  the  four  estimation  routines,  a 
regression  was  then  performed  of  y  on  X  and  residuals  calculated  using, 


,  J*J' m ijij ij  'ir  ;  m 


'L'  ns1*  i  'jr.’f :  f.v;  ■■.  ■■ 


TABLE  I 

ESTIMATES  OF  AUTOCORRELATION  COEFFICIENT 


p 

DW 

TN 

PW 

BM 

.2 

.19 

.19 

.19 

.19 

.6 

.59 

.60 

.58 

.60 

.8 

.78 

.80 

.77 

.80 

The  values  of  the  residuals  were  then  sent  to  each  estimation  routine.  Estimates  of 
p  and  p  were  determined  for  values  of  p  equal  to  .2,  .3,  .4,  .5,  .6,  .7,  .8,  and  .9.  Each 
estimate  was  replicated  1000  times  for  the  sample  sizes  of  20  and  50. 


F.  MEASURES  OF  EFFECTIVENESS 

In  order  to  compare  the  performances  of  the  estimators,  two  MOE’s  were  used. 

A 

The  mean  square  error  (MSI:)  of  p  was  estimated  for  each  estimator.  This  represents 
the  expected  squared  error  made  in  estimating  p.  The  following  computational  form  of 
MSE  was  used, 


y  (p  -  p,)2/iooo 


i=  1 . 1000. 


(cqn  3.3) 


The  successive  values  of  MSE  of  p  were  then  plotted  against  the  actual  p  to 
examine  performance  over  the  range  of  p. 


The  second  MOE  examined  the  relative^clTicicncies  of  the  regression  coefficient  a^ 
defined  in  (Ref.  7:p.  7|.  A  ratio  of  MSE  of  P  for  a  particular  estimate  to  the  MSE  of  P 


IV.  RESULTS  AND  CONCLUSIONS 

A.  GENERAL 

The  major  emphasis  of  this  thesis  was  to  examine  the  performance  of  four 
estimators  of  the  autocorrelation  coefficient,  p,  for  a  first  order  autoregressive  process. 
The  estimators  examined  were  Durbin-Watson,  Theil-Nagar,  Prais-Winsten,  and 
Beach-MacKinnon. 

A  Monte-Carlo  simulation  was  performed  for  the  following  values  of  p:  .2,  .3,  .4, 
.5,  .6,  .7,  .8,  and  .9.  Each  run  was  replicated  1000  times  for  sample  sizes  of  20  and  50. 
The  results  are  recorded  in  Table  II  .  Irrespective  of  sample  size,  each  of  the  methods 
underestimate  the  true  value  of  p  but  as  the  number  of  observations  is  increased  from 
20  to  50  the  bias  reduces.  As  was  expected,  no  one  estimator  uniformly  outperforms 
the  others.  In  both  sample  sizes,  the  two  stage  estimators  (Durbin-Watson  and 
Theil-Nagar)  achieve  better  results  for  small  p.  As  the  value  of  p  increases,  the 
iterative  methods  (Prais-Winsten  and  Beach-MacKinnon)  perform  best.  With  T=20 
this  transition  occurs  at  p  =  .6  while  at  50  observations  it  occurs  earlier  at  p  =  .4. 

The  discussion  of  the  results  will  be  divided  into  two  sections.  The  measures  of 
effectiveness,  as  defined  in  Chapter  3  will  first  be  applied  to  the  simulation  results  for 
T=20.  This  will  be  followed  by  an  identical  approach  when  the  sample  size  is 
increased  to  50. 

I 

B.  SAMPLE  SIZE  20 

Since  performance  of  an  estimator  is  roughly  indicated  by  its  mean  and  variance, 
mean  square  error  (MSE)  of  each  p1  over  the  entire  sample  size  was  estimated.  The 
results  of  these  calculations  arc  presented  in  Table  III  along  with  a  plot  of  MSE  of  p'' 
versus  actual  values  of  p  in  Figure  4.1  .  They  again  indicate  that  the  Theil-Nagar  and 
Durbin-Watson  estimators  arc  better  for  smaller  values  of  p  (p<.6)  and  as  p  increases 
the  Prais-Winsten  p  emerges  as  the  best.  On  the  basis  of  Figure  4.1  alone, 
Beach-MacKinnon's  performance  is  clearly  inferior.  However,  in  examining  the 
efficiency  of  each  estimator  in  Table  IV,  Beach-MacKinnon  proves  to  be  the  most 
efficient  in  estimating  P  over  the  widest  range  of  p.  The  tie  in  Figure  4.1  between 
Theil-Nagar  and  Durbin-Watson  is  resolved  in  Table  IV  with  Theil-Nagar  proving  to 
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TABLE  II 

ESTIMATES  OF  AUTOCORRELATION  COEFFICIENT 


Sample  Size  20 


p 

DW 

TN 

PVV 

BM 

.2 

.158 

.162 

.113 

.107 

.3 

.234 

.239 

.205 

.193 

.4 

.310 

.316 

.296 

.279 

.5 

.385 

.392 

.387 

.365 

.6 

.460 

.467 

.478 

.450 

.7 

.533 

.542 

.567 

.535 

.8 

.603 

.613 

.655 

.617 

.9 

.667 

.681 

.741 

.697 

Sample  Size  50 

P 

DW 

TN 

pw 

BM 

.2 

.173 

.161 

.170 

.178 

.3 

.279 

.253 

.270 

.270 

.4 

.360 

.340 

.360 

.359 

.5 

.451 

.432 

.463 

.453 

.6 

.544 

.523 

.559 

.547 

.7 

.630 

.610 

.651 

.640 

.8 

.726 

.700 

.747 

.731 

.9 

.812 

.796 

.839 

.820 

be  uniformly  more  efficient  than  Durbin-Watson.  Table  IV  also  demonstrates  that  for 
p  =  .2  OLS  is  at  least  as  efficient  as  three  of  the  four  estimators. 

C.  SAMPLE  SIZE  50 

The  results  of  the  \1SE  calculations  for  T=50  arc  recorded  in  Table  V  along 
with  a  plot  of  MSE  of  p  versus  the  actual  values  of  p  in  Figure  4.2.  The 
Durbin-Watson  and  Theil-Nagar  estimators  again  perform  the  best  for  smaller  values 
of  p  (p  <  .4)  and  as  p  increases  the  Bcach-McKinnon  and  Prais-Winsten  estimators  of 
p  contain  the  smallest  MSE. 

Once  again  even  though  the  Prais-Winsten  p  has  a  smaller  MSE  than 
Beach-MacKinnon,  Table  VI  illustrates  that  Beach-McKinnon  is  a  uniformly  more 
efficient  estimator  of  the  slope  coefficient,  for  the  smaller  values  of  p  (p<.4) 
Theil-Nagar  is  more  efficient  than  Durbin-Watson.  Table  VI  also  illustrates  that  OLS 
is  at  least  as  efficient  as  any  of  the  other  estimators  when  p  is  small. 

D.  SUMMARY 

As  was  found  in  previous  studies  when  autocorrelation  is  present  only  to  a  slight 
degree  (p<.2)  the  OLS  estimator  provides  an  efficient  estimate  for  the  regression 
coefficient,  p.  As  the  process  becomes  more  significant  however,  all  the  estimators 
outperform  the  OLS  solution.  In  both  sample  sizes  the  performance  of  Theil-Nagar 
and  Durbin-Watson  are  nearly  identical  witli  respect  to  the  MSE  of  p.  However,  when 
efficiency  of  the  slope  coefficient  estimate  is  examined,  Theil-Nagar  proves  to  be  the 
better  2  stage  estimator.  Park  and  Mitchell  [Ref.  7:p.  4]  found  that  Prais-Winsten 
performs  better  in  estimating  p.  The  results  presented  here  tend  to  dispute  that 
finding.  for  while  Prais-Winsten  has  a  uniformly  smaller  MSE  of  p, 
Beach-MacKinnon  provides  the  most  efficient  estimator  of  P  .  Spitzcr,  on  the  other 
hand  [Ref.  ll:p.  44),  which  ranked  two  stage  estimators  as  being  the  best  for  values  of 
p  between  .2  and  .5,  mirrors  the  results  produced  here.  Apriori  knowledge  of  the 
neighborhood  of  p  will  be  helpful  in  selecting  the  appropriate  estimation  method.  For 
both  sample  sizes  Theil-Nagar  appears  to  be  the  best  for  small  values  of  p. 
Beach-MacKinnon,  while  containing  a  larger  bias  for  p  than  docs  Prais-Winsten,  is  a 
much  more  efficient  estimator  of  the  slope  coefficient  for  larger  values  of  p. 
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TABU;  IV 

F.ITTCIHXCY  OF  lUiGRliSSIOX  COIITlCHiXT  IlSTIMA TIiS 


Sample  Si/e  20 


l> 

MSFp  (I)\V ) 

MSIip  (TN) 

MSI  ip  (l’\V) 

M Slip  (BM) 

M Slip  (OI.S) 

MSIiP  (OI.S) 

MSIip  (OI.S) 

MSIip  (OI.S) 

2 

1.004 

.9794 

1.035 

1.041 

3 

.9228 

.8967 

.9442 

.9515 

4 

.8218 

.7929 

.8325 

.8342 

5 

.7082 

.6751 

.7024 

.6959 

6 

.5864 

.5484 

.5652 

.5515 

7 

.4610 

.4207 

.4329 

.4135 

8 

.3359 

.3020 

.3093 

.2870 

<) 

.2253 

.2087 

.2077 

.1892 

TAB I-II  VI 

HFFK'lIiNCY  OF  RI:G  RUSS  ION  COUI  I  ICIUNT  USTIMATUS 
Sample  Size  50 


p 

MSl.p  (D\V) 

MSl.p  (TN) 

MS  Up  (P\V) 

MSCp  (BM) 

MSnp  (OLS) 

MSI-p  (OI.S) 

Msr-p  (oi.s) 

MSFp  (OI.S) 

.2 

1.073 

1.041 

1.046 

1.058 

.3 

.9985 

.94S2 

.9714 

.9562 

.4 

.8850 

.8255 

.8635 

.8255 

.5 

.7452 

.6859 

.7282 

.6825 

.6 

.5920 

.5420 

.5870 

.5406 

.7 

.4366 

.4020 

.4453 

.4020 

.8 

.2889 

.2690 

.3067 

.2700 

.9 

.1589 

.1505 

.1738 

.1505 

2' 


APPENDIX 
PROGRAM  LISTINGS 

This  appendix  contains  listings  of  the  programs  utilized  in  the  analysis  performed 
herein.  All  of  the  functions  are  written  in  FORTRAN  and  contain  the  necessary 
documentation.  The  Monte  Carlo  simulation  was  performed  using  the  Advanced 
Simulation  and  Statistics  Package  [Ref.  12]  developed  by  P.  A.  Lewis.  Since  the 
package  only  allows  for  the  simultaneous  comparision  of  3  estimators,  2  functions  were 
developed  for  each  sample  size.  The  first,  SIMS  generates  estimates  for 
Durbin- Watson,  Theil-Nagar,  and  Prais-Winsten  for  a  sample  size  of  20.  SI  MSA 
meanwhile,  generates  estimates  for  Bcach-MacKinnon  for  the  identical  sample  size. 
Routines  for  Durbin-Watson  and  Theil-Nagar  were  included  in  SI  MSA  to  ensure  the 
results  were  comparable  to  SIMS.  SIMSB  and  SIMSC  perform  in  a  similar  fashion  for 
sample  size  of  50  and  therefore  were  not  included.  The  Advanced  Simulation  and 
Statistics  Package  computes  the  mean  square  error  of  'p'  for  each  estimator 
automatically.  The  mean  square  error  for  the  ji  estimates  was  accomplished  by  the 
MSEB  runction. 


SIMS 

DIMENSION  EHAT(20) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 

COMMON  /DATA1/  IX1A.RH0 

REAL* 4  Y( 5000 ) , YMIN , YMAX , PMEAN( 3 ) 

CHARACTER*80  TI,T2,T3 

INTEGER  N, M,NE(8) ,L,D,RG,SEI ,SVS, NEST, NSR.IXl, 1X2, 1X3 
EXTERNAL  DATGEN,  DURWAT,  BEAMAC,  PRAWIN,  LSEB,  DCALC,  TRANSF 
EXTERNAL  LN0RM , S IMTBD , GMPRD 
NR=20 
T=20 
K=2 
C 
C 
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OPEN  (UNIT=19,FIL£='MONICAl ) 

C  OPEN  (UNIT=21 , FILE=' MARGE' ) 

C  OPEN  (UNIT=51,FILE='AMBROSE') 

C  OPEN  (UNIT=41,FILE='DAT2') 

OPEN  ( UN  I T=6 1 , FI LE= 1  DAT 3 1 ) 

READ  (19,*)  ANS 

10  READ( 19 , * , END=999)  N,M, L,D,RG,SEI , SVS ,NEST,NSR 

READ( 19 , * )YMIN , YMAX 
READ( 19 ,*)  (NE( I) , 1=1 , L) 

READ( 19 , 120)  IX1.IX2.IX3 
120  FORMAT( 15 , IX, 15, IX, 15) 

READ  (19,115)  T1 
115  F0RMAT(A80) 

READ( 19,115)  T2 
READ  ( 19 , 1 15 )T3 
READ(19,*)  (PMEAN(I) ,1=1,3) 

READ(19,*)  RHO 
READ(19,61)IX1A 
61  F0RMAT(I5) 

C 

C  CALL  FOR  SIMTBD 

C 

CALL  SIMTBD  (IX1,IX2,IX3,Y,N,M,NE,L,D,NSR,RG,SEI,SVS, 

*YMIN , YMAX , NEST , DATGEN , DURWAT , T 1 , DATGEN , BEAMAC , T2 , DATGEN , PRAWI N , T3 , 
*PMEAN) 

GO  TO  10 

999  WRI TE( 6 , * ) 1  END  OF  DATA  INPUT' 

STOP 
END 

Q  I**************************************************************** 

C  *************  OATA  GENERATION  SUBROUTINE  ********************** 

p  *******************rt*Ar*****rt**rtrt**rt**rtrt*rt**********rt**rtrt*******rt* 


SUBROUTINE  DATGEN  ( 1X1 ,EHAT,NR) 

DIMENSION  BHAT(2) ,YSTAR(20) ,R2(20)  ,U(20) , 


*E(20) ,YHAT(20) ,EHAT(20) ,XSTAR(20,2) 

*,Y1(20),X(20,2),V(20) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 
COMMON  /DATA1/  IX1A.RHO 
INTEGER  I X 1 , 1 X 1A , NR 
C 
C 

C  GENERATE  THE  RANDOM  ERROR 

C 

CALL  SNOR  (IX1,U,NR,1,0) 

C 

C  ADJUST  THE  VARIANCE  OF  R.  E.  IAW  BEACH  AND  MACKINNON( 1978) 

DO  38  1=1 ,T 

V(I)=U(I)\  06 
38  CONTINUE 

C 

C  GENERATE  THE  ERROR  FOR  THE  STAND  LINEAR  MODEL 

C 

E(1)=V(1)/(1-(RH0**2))**0. 5 
DO  31  J=2 ,T 

E( J)=RHQ*E( J-l)+V( J) 

31  CONTINUE 
C 

C 

C  GENERATE  THE  EXPLANATORY  VARIABLES  IAW  RAO  AND  GRILITCHES  (1969) 

C 

DO  32  1=1,20 
X(I,1)=1 

32  CONTINUE 

C  CHANGE  1X1  IN  ORDER  TO  AVIOD  COLLINEARITY 

C  IX1A=IX1+19 

CALL  SNOR( IX1A.R2 ,NR, 1 ,0) 

DO  33  J=1 ,20 

X(J  ,2)=R2(J)*.  25 

33  CONTINUE 
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THE  TRUE  BETA  EQUALS  1,1 


C 
C 
C 

C  GENERATE  THE  INDEPENDENT  VARIABLE 
C 

DO  35  1=1,20 

Y 1 ( I  )=( X ( 1 , 1)+X( 1 ,2))+E( I) 

35  CONTINUE 

C 

C  GENERATE  THE  LEAST  SQUARES  ESTIMATOR  FOR  BETA 

C 

CALL  LSEB(X,Y1,BHAT) 

C  PRINT  LSEB  TO  A  FILE 

IF(ANS  .EQ.  2)  WRITE(61 ,201)BHAT 
201  F0RMAT(F11.  8,2X,F11. 8) 

C 

C  GENERATE  YHAT 

C 

DO  100  1=1,20 

YHAT(I)=X(I , 1)*BHAT( 1)+X( I ,2)*BHAT(2) 

100  CONTINUE 

C 

C 

C  GENERATE  EHAT 

C 

DO  50  1=1,20 

EHAT ( I )=Y 1 ( I ) -YHAT ( I ) 

50  CONTINUE 

C 

C 

RETURN 

END 

C  *********************  0URBIN  WATSON  ******* 
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C  THIS  FUNCTION  COMPUTES  THE  DURBIN-WATSON  ESTIMATE  OF  RHO 
REAL  FUNCTION  DURWAT  (EHAT,NR,WI) 

DIMENSION  EHAT(20) ,X(20,2),Y1(20) ,XSTAR1(20,2) ,YSTAR1(20) , BHAT 1(2) 
COMMON  /MYDATA/  K,T,ANS,Y1,X 
CALL  DCALC  (EHAT,T,D) 

DURWAT=l-D/2 

C 

CALL  TRAN$F(X,Y1 .DURWAT, XSTAR1 ,YSTAR1) 

CALL  LSEB  (XSTAR1 , YSTAR1 ,BHAT1) 

IF  (ANS  .EQ.  1  )  WRITE(21 ,701)  BHAT1 
701  F0RMAT(F11.8,2X,F11.8) 

C 

C 

C 

RETURN 

END 

C 

C 

C  ********************  THEIL  NAGAR  ****************************** 

C  THIS  FUNCTION  COMPUTES  THE  THEIL-NAGAR  ESTIMATE  OF  RHO 
C  REAL  FUNCTION  THENAG  (EHAT.NR.WI) 

DIMENSION  EHAT(20),YSTAR2(20),XSTAR2(20,2),BHAT2(2) 

*  ,Y1(20) ,X(20 ,2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 
CALL  DCALC  (EHAT,T,D) 

THENAG=( (T**2)*( l-D/2)+K**2)/(T**2-K**2) 

CALL  TRANSF( X , Y 1 , THENAG , XSTAR2 , YSTAR2 ) 

CALL  LSEB  ( XSTAR2 , YSTAR2 ,BHAT2) 

IF  (ANS  .EQ.  1  )  WRI TE( 3 1 ,801)  BHAT2 
801  F0RMAT(F11.  8,2X,F11.  8) 

RETURN 

END 

C  ***********************  PRAIS  WINSTEN  ************************* 

C  THIS  FUNCTION  COMPUTES  THE  PRAIS-WINSTEN  ESTIMATE  OF  RHO 
REAL  FUNCTION  PRAWIN( EHAT.NR.WI) 


DIMENSION  EHAT3(2Q) ,YHAT3(20) , YSTAR3(20) ,BHAT3(2) , 
*EHAT(20) , XSTAR3( 20,2) 

^ , Y1 ( 20 ) ,X(20,2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 
N=0 

RH03=0 

N=N+1 

CALL  TRANS F  (X,Y1,RH03,XSTAR3,YSTAR3) 

CALL  LSEB  (XSTAR3 ,YSTAR3 ,BHAT3) 

GENERATE  YHAT3 
DO  83  1=1,20 

YHAT3( I)=X( I , 1)*BHAT3(1)+X( I ,2)*BHAT3(2) 
CONTINUE 
DO  4  1=1, T 

EH AT  3(I)=Y1(I) — YH AT  3(1) 

CONTINUE 

RH0NUM=0 
RHODEN=0 
DO  5  1=2, T 

RHONUM=RHQNUM+( EHAT3( I)*EHAT3( 1-1 ) ) 

CONTINUE 


DO  6  1=2, T-l 

RHODEN=RHODEN+( EHAT3( I)**2) 

CONTINUE 

PRAWIN=RHONUM/RHODEN 

CHECK  FOR  PRAWIN  WHICH  ARE  OUT  OF  BOUNDS 
IF(PRAWIN.  GE.  1)THEN 
PRAWIN=0. 99999 
ELSE  IF  (PRAWIN. LE. -1)THEN 
PRAWIN=-0.  99999 
END  IF 

COMPARISION  OF  RH03  AND  PRAWIN  IF  DIFF  . LT.  0.0001  THEN  END 
IF(ABS(RH03-PRAWIN).  GT. .  0001)THEN 
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RH03=PRAWIN 
GO  TO  98 
ELSE 

PRAWIN=PRAWIN 
END  IF 

C  IF  (ANS  .EQ.  1  )  WRITE(41 ,901)  BHAT3 
CO 1  F0RMAT(F11.  8,2X,F11.  8) 

RETURN 

END 

C 

C  THE  FOLLOWING  SUBROUTINES  AID  IN  THE  COMPUTATION  OF  THE  FOUR 
C  ESTIMATORS  OF  RHO. 

C  *************  SUBROUTINE  LSEB  ******************************* 

C  SUBROUTINE  LSEB  WILL  COMPUTE  THE  LSE  OF  B 
C 

SUBROUTINE  LSEB( X , Y1 , BHAT) 

DIMENSION  BHAT(2) ,Y1(20),X(20,2) ,XTRNSP(2 ,20) ,XI(2,2),H(2,20), 
*XPRIX(2,2) 

C  X  TRANSPOSE 

DO  40  1=1,20 
DO  41  J=1 ,2 

XTRNSP(J , I )=X( I , J ) 

41  CONTINUE 

40  CONTINUE 

C  MULTIPLY  X  TRANSPOSE  AND  X 

CALL  GMPRD(XTRNSP,X,XPRIX,2,20,2) 

C  CALCULATE  INVERSE  OF  X  PRIME  X 

DETR=1/(XPRIX(1,1)*XPRIX(2,2)-XPRIX(1,2)*XPRIX(2,1)) 

X I  ( 1 , 1)=DETR*XPRIX(2 ,2) 

XI ( 1 ,2)=DETR*(-XPRIX(  1,2)) 

XI(2 , 1)=DETR*(-XPRIX(2, 1)) 

XI(2 , 2)=DETR*XPRIX(  1,1) 

C  MULTIPLY  INVERSE  AND  TRANSPOSE 

CALL  GMPRD(XI , XTRNSP ,H ,2 ,2 ,20) 

DO  99  1=1,2 


BHAT ( I )=H( I,1)*Y1(1)+H(I,2)*Y1(2)+H(I,3)*Y1(3) 

*+H(  I,4)*Y1(4)+H(I,5)*Y1(5) 

*+H(1 ,6)*Y1(6)+H(I ,7)ytYl(7 )+H( 1 ,8)*Y1(8)+H( 1 ,9)*Y1(9) 

*+H( I , 10)*Y1(10)+ 

*H(I , 11)*Y1(11)+H( I , 12)*Y1( 12)+H( I ,13)*Y1( 13)+H( I , 14)*Y1( 14) 
*+H(I.15)*Yl(15)+ 

*H(1 ,16)*Y1( 16)+H( 1 ,17)*Y1( 17)+H( 1 ,18)*Y1( 18)+H( I,19)*Y1(19) 
*H( I ,20)*Y1(20) 

CONTINUE 

RETURN 

END 

*************  SUBROUTINE  DCALC  *************************** 
SUBROUTINE  DCALC  WILL  COMPUTE  THE  DURBIN  STATISTIC  D 

SUBROUTINE  DCALC(EHAT,T,D) 

DIMENSION  Dl(20) ,D2(20) ,EHAT(20) 

DNUM=0 
DDEN=0 
DO  1  1=2, T 

D 1 ( I — 1 )=( EHAT ( I ) — EHAT ( I-l))**2 
DNUM=DNUM+D1( 1-1)  - 

CONTINUE 
DO  2  J=1 ,T 

D2( J)=EHAT( J)**2 
DDEN=DDEN+D2( J) 

CONTINUE 

D=DNUM/DDEN 

RETURN- 

END 

*****************  SUBROUTINE  TRANSF  ********************** 

SUBROUTINE  TRANSF  IS  DESIGNED  TO  TRANSFORM  THE  X'S  AND  Y'S 
ACCORDING  TO  THE  LEAST  SQUARES  RULE. 

SUBROUTINE  TRANSF(X,Yl,RHOHAT,XSTAR,YSTAR) 


DIMENSION  Y 1 ( 20 ) ,YSTAR(20) ,X(2Q,2) ,X$TAR(20,2) 

K=2 

T=20 

Y  TRANSFORM 

YSTAR( 1)=((1-(RH0HAT**2))**0.  5)*Y1(1) 

DO  7  1=2,20 

YSTAR( I )=Y1( I )-( RH0HAT*Y1( 1-1 ) ) 

CONTINUE 
X  TRANSFORM 

XSTAR( 1 , 1)=( I-(RH0HAT**2))**0.  5 
DO  9  J=2 , K 

XSTAR(1 , J)=(( 1-(RHOHAT**2))**0. 5)*X(1,J) 
CONTINUE 
DO  11  L=2,T 

XSTAR( L, 1 )=l-RHOHAT 
CONTINUE 
DO  12  1=2, T 
DO  13  J=2,K 

XSTAR(I , J)=X( I , J)-RHQHAT*X( 1-1 , J) 
CONTINUE 
CONTINUE 
RETURN 


SIMSA 

THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  RUN  COMPUTE  THE  FOLLOWING 
ESTIMATORS  (DW  TN  BM)  FOR  A  SAMPLE  SIZE  OF  20 
DIMENSION  EHAT(20) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 
COMMON  /DATA1/  IX1A.RHO 
REALM  Y( 5000 ) , YMIN , YMAX , PMEAN( 3 ) 

CHARACTER*80  T1.T2.T3 

INTEGER  N, M,NE(8),L,D,RG,SEI,SVS, NEST, NSR.IXl, 1X2, 1X3 
EXTERNAL  DATGEN,  DURWAT,  THENAG,  BEAMAC,  LSEB,  DCALC ,  TRANS F 
EXTERNAL  LNORM,$IMTBD,GMPRQ 
NR=20 


OPEN  ( UNIT=19 , FI LE= ' MONICA ' ) 

OPEN  (UNIT=51 , FI LE=‘ AMBROSE' ) 

READ  (19,*)  ANS 

READ( 19 , * , END=999)  N,M,L,D,RG,SEI ,SVS,NE$T,NSR 
READ( 19,*) YMIN, YMAX 
READ( 19 ,*)  (NE( I ) , 1=1 , L) 

WRITE  (22,105)  (NE( I) , 1=1 , L) 

F0RMAT(8I4) 

READ( 19 , 120)  IX1.IX2.IX3 
F0RMAT( 15, IX , 15 , IX , 15) 

READ  (19,115)  T1 
F0RMAT(A80) 

READ( 19,115)  T2 
READ  ( 19 , 115)T3 
READ( 19,*)  ( PMEAN( I ) , 1=1 ,3) 

READ( 19 ,*)  RHO 
READ( 19 ,61)IX1A 
FORMAT(  15) 


c 

C  CALL  FOR  SIMTBD 
C 

CALL  SIMTBD  ( 1X1 , 1X2, 1X3 , V.N.M.NE, L.D.NSR, RG.SEI ,SVS , 

*YMIN, YMAX , NEST , DATGEN , DURWAT ,T1 , DATGEN , THENAG , T2 , DATGEN , BEAMAC , T3 , 
*PMEAN) 

GO  TO  10 

999  WRITE(6 ,*) ' END  OF  DATA  INPUT' 

STOP 

END 

Q  *************************  ******  ********************************* 

C  *************  DATA  GENERATION  SUBROUTINE  ********************** 

Q  ***************************************************************** 

SUBROUTINE  DATGEN  ( I X 1 , EHAT , NR) 

DIMENSION  BHAT(2) ,YSTAR(20) ,R2(20) , U( 20 ) , 

*E(20) ,YHAT(20) ,EHAT(20) ,XSTAR(20,2) 

*,Y1(20) ,X(20,2) 

COMMON  /MYDATA/  K.T.ANS.Yl.X 
COMMON  /DATA1/  IX1A.RHO 
INTEGER  1X1 , IX1A.NR 
C 
C 

C  GENERATE  THE  RANDOM  ERROR 

C 

CALL  SNOR  (1X1 ,U,NR, 1 ,0) 

C 

C 

C  GENERATE  THE  ERROR  FOR  THE  STAND  LINEAR  MODEL 

C 

E( 1)=U( l)/( l-(RHO**2) )**0. 5 
DO  31  J=2,20 

E(U)=RHO*E(J-l)+U( J) 

31  CONTINUE 

C 

C 
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GENERATE  THE  EXPLANATORY  VARIABLES  IAW  RAO  AND  GRI LITCHES  (1969) 

DO  32  1=1,20 
X(I,1)=1 
CONTINUE 

CHANGE  1X1  IN  ORDER  TO  AVIOD  CQLLINEARITY 
IX1A=IX1+19 

CALL  SNOR( IX 1A , R2 , NR , 1 , 0 ) 

DO  33  J=1 ,20 

X(J,2)=R2(J)*.  25 
CONTINUE 


THE  TRUE  BETA  EQUALS  1,1 


GENERATE  THE  INDEPENDENT  VARIABLE 
DO  35  1=1,20 

Yl( I )=(X( 1 , 1)+X( 1 ,2))+E(  I) 
CONTINUE 


GENERATE  YHAT 

CALL  LSEB(X ,Y1 ,BHAT) 

BHAT( 1)=1.  3 
BHAT(2)=1.  1 
DO  100  1=1,20 

YHAT( I )=X( I , 1)*BHAT( 1)+X( I ,2)*BHAT(2) 
CONTINUE 


GENERATE  EHAT 
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C 


C 

c 


c 

c 

c 


c 

c 

c 

c 


c 

c 


DO  50  1=1,20 

EHAT( I)=Y1(I) -YHAT ( I ) 
CONTINUE 


RETURN 

END 

*********************  DURBIN  WATSON  *************************** 
REAL  FUNCTION  DURWAT  (EHAT,NR,WI) 

DIMENSION  EHAT(20) ,X(20,2),Y1(20) ,XSTAR1(20,2) , YSTAR1 ( 20 ) , BHAT 1(2) 
COMMON  /MYDATA/  K,T,AN$,Y1,X 
CALL  DCALC  (EHAT,T,D) 

DURWAT=l~D/2 

CALL  TRANS F( X ,Y1 , DURWAT, XSTAR1 ,YSTAR1) 

CALL  LSEB  (XSTAR1 ,YSTAR1 ,BHAT1) 


RETURN 

END 


********************  THE I L  NAGAR  ****************************** 


REAL  FUNCTION  THENAG  (EHAT,NR,WI) 

DIMENSION  EHAT(20) , YSTAR2( 20 ) , X STAR2 ( 20 ,2) ,BHAT2(2) 

*,Y1(20),X(20,2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 
CALL  DCALC  ( EHAT.T.D) 

THENAG=((T**2)*( l-D/2)+K**2)/(T**2-K**2) 

RETURN 
END  ' 

***********************  BEACH  MACKINNON  ************************* 
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REAL  FUNCTION  BEAMAC( EHAT , NR , WI ) 

DIMENSION  EHAT4(20) ,YHAT4(20) ,YSTAR4(20) ,BHAT4(2) , 
*Y1(20) ,EHAT(20) ,X(20,2) ,XSTAR4(20 ,2) 

COMMON  /MYDATA/  K,T,ANS,Y1,X 
N=0 

RH04=0 

N=N+1 

CALL  TRANSF  (X ,Y1 ,RH04,XSTAR4 ,YSTAR4) 

CALL  LSEB  (XSTAR4,YSTAR4,BHAT4)  . 

BHAT4(1)=1.  0 
BHAT4(2)=1.  0 
GENERATE  YHAT4 
DO  83  1=1,20 

YHAT4(I)=X(I,1)*BHAT4(1)+X(I,2)*BHAT4(2) 
CONTINUE 
DO  4  1=1, T 

EHAT  4(1 )=Y 1(1 ) — YHAT  4(1) 

CONTINUE 
SUM3=Q 
SUM2=0 
SUM1=0 
DO  71  1=2, T 

SUM1=SUM1+(EHAT4( I )*EHAT4( 1-1) ) 

CONTINUE 

DO  72  1=2, T 

SUM2=$UM2+(EHAT4( I-l)**2) 

CONTINUE 

DO  73  1=2, T 

SUM3=SUM3+(EHAT4(I)**2) 

CONTINUE 

DEN0M=(T-1)*(SUM2-(EHAT4( 1)**2) ) 


A=(-(T-2)*SUM1)/DEN0M 

c 

8=(((T“1)*(EHAT4( 1)**2))-(T*SUM2)-SUM3)/DEN0M 
C 

C=(T*SUM1)/DEN0M 

C 

SMALP=B-((A**2)/3) 

C 

SMALQ=C-((A*B)/3)+((2*(A**3))/27) 

C 

THETA=AC0S((SMALQ*(27**.  5))/(2*SMALP*((-SMALP)**0.  5))) 

C 

C  BEAMAC  IS  THE  ITERATIVE  RHO  FOR  THIS  PROCEEDURE 

BEAMAC=(-2*((-SMALP/3)**0.  5))*COS((THETA/3)+(3. 1412/3) )-( A/3) 
C  CHECK  FOR  BEAMAC  WHICH  ARE  OUT  OF  BOUNDS 
I F( BEAMAC.  GE.  1)THEN 
BEAMAC=0.  99999 
ELSE  IF  (BEAMAC. LE.-1)THEN 
BEAMAC=-0.  99999 
END  IF 

C  COMPARISION  OF  RH04  AND  BEAMAC  IF  DIFF  . LT.  O.OOOl  THEN  END 
IF(ABS(RH04-BEAMAC).  GT. .  0001)THEN 
RH04=BEAMAC 
GO  TO  98 
ELSE 

BEAMAC=BEAMAC 
END  IF 

IF  (ANS  .EQ.  2)  WRITE  (51,901)  BEAMAC 
901  FORMAT( F15. 11) 

RETURN 

END 

C 

C  THE  FOLLOWING  SUBROUTINES  AID  IN  THE  COMPUTATION  OF  THE  FOUR 
C  ESTIMATORS  OF  RHO. 


*************  SUBROUTINE  LSEB  ******************************* 
SUBROUTINE  LSEB  WILL  COMPUTE  THE  LSE  OF  B 

SUBROUTINE  LSEB(X ,Y1 ,BHAT) 

DIMENSION  BHAT(2) , Yl(20) ,X(20 ,2) ,XTRNSP(2,20) ,XI(2,2) ,H(2 ,20) , 
*XPRIX(2,2) 

X  TRANSPOSE 
DO  40  1=1,20 
DO  41  J=1 ,2 

XTRNSP( J , I )=X( 

CONTINUE 

CONTINUE 

MULTIPLY  X  TRANSPOSE  AND  X 
CALL  GMPRD(XTRNSP,X,XPRIX,2,20,2) 

CALCULATE  INVERSE  OF  X  PRIME  X 
DETR=1/(XPRIX(1,1)*XPRIX(2,2)-XPRIX(1,2)*XPRIX(2,1)) 

XI ( 1 , 1)=DETR*XPRIX(2,2) 

XI( 1 ,2)=DETR*(-XPRIX( 1,2)) 

XI(2,1)=DETR*(-XPRIX(2,1)) 

XI(2,2)=DETR*XPRIX( 1,1) 

MULTIPLY  INVERSE  AND  TRANSPOSE 
CALL  GMPRD(XI,XTRNSP,H, 2,2,20) 

DO  99  1=1,2 

BHAT(I)=H(I,1)*Y1(1)+H(I,2)*Y1(2)+H(I,3)*Y1(3) 

*+H( I,4)*Y1(4)+H(I,5)*Y1(5) 

*+H( I,6)*Y1(6)+H(I,7)*Y1(7)+H(I,8)*Y1(8)+H(I,9)*Y1(9) 
*+H(I,10)*Yl(10)+ 

*H( I , 11)*Y1(11)+H(I , 12)*Y1( 12)+H( I , 13)*Y1( 13)+H( I,14)*Y1(14) 
*+H(I,15)*Yl(15)+ 

*H(I,16)*Y1(16)+H(I, 17)*Y1( 17)+H( I,18)*Y1(18)+H(I,19)*Y1(19)+ 
*H(I ,20)*Y1(20) 

CONTINUE 

RETURN 

END 


************* 


SUBROUTINE  DCALC 


***************************** 


C  SUBROUTINE  DCALC  WILL  COMPUTE  THE  DURBIN  STATISTIC  D 
C 

SUBROUTINE  DCALC (EHAT,T,D) 

DIMENSION  Dl(20) ,D2(20) ,EHAT(20) 

DNUM=0 
DDEN=0 
DO  I  1=2, T 

D1(I-1)=(EHAT(I)-EHAT(I-1))**2 
DNUM=DNUM+D1( 1-1) 

1  CONTINUE 

DO  2  J=1,T 

D2(J)=EHAT(J)**2 

DDEN=DD£N+D2(J) 

2  CONTINUE 

D=DNUM/DDEN 

RETURN 

END 

C 

C  *****************  SUBROUTINE  TRANSF  ********************** 

C 

C  SUBROUTINE  TRANSF  IS  DESIGNED  TO  TRANSFORM  THE  X'S  AND  Y'S 
C  ACCORDING  TO  THE  LEAST  SQUARES  RULE. 

SUBROUTINE  TRANSF(X ,Y1 ,RHOHAT,XSTAR,YSTAR) 

DIMENSION  Y I ( 20 ) ,YSTAR(20) ,X(20 ,2) ,XSTAR(20 ,2) 

K=2 

T=20 

C  Y  TRANSFORM 

YSTAR(1)=((1-(RH0HAT**2))**0.  5)*Y1(1) 

DO  7  1=2,20 

YSTAR( I )=Y1( I )-( RH0HAT*Y1( 1-1) ) 

7  CONTINUE 

C  X  TRANSFORM 

XSTAR(  1 , 1)=(  1-(RH0HAT**2)  )*'*0.  5 
DO  9  J=2,K 

XSTAR( 1 , J )=( ( 1-(RH0HAT**2) )**0. 5)*X(1,J) 
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CONTINUE 
DO  11  L=2,T 

XSTAR(L,l)=l-RHOHAT 
CONTINUE 
DO  12  1=2, T 
DO  13  J=2,K 

XSTAR( I , J)=X( I , J)-RHOHAT*X( 1-1 ,J) 
CONTINUE 
CONTINUE 
RETURN 


THIS  PROGRAM  IS  DESIGNED  TO  CALCULATE  THE  MEAN  SQUARE  ERROR  OF 
THE  BETA  VECTOR 

DIMENSION  Bl( 5000) ,62(5000) , B3( 5000 ) , B4( 5000 ) ,B5(5000) , B6 ( 5000 ) , 
*B7(5000) ,88(5000), B9( 5000) ,B10(5000) ,BX(5000) ,BY(5000) 

OPEN  (UNIT=21 , FI LE= ' DAT1 1 ) 

OPEN  (UNIT=31 , FI LE= 1  DAT 2 ' ) 

OPEN  ( UN IT=4 1 , FILE='DAT3 ' ) 

OPEN  (UNIT=51 , FI LE= 1  DAT 4 ' ) 

OPEN  (UNIT=61,FILE='DAT5') 

COUNT=1000 

READ(21 , 900 ) ( B1 ( I ) ,B2( I ) ,  1=1,1000) 

CALL  MSEBET  (B1.B2, COUNT, XMSEDW) 

READ(31 ,900)(B3( I ) ,84( I) ,  1=1,1000) 

CALL  MSEBET  (B3,B4,C0UNT,XMSETN) 

READ(41 ,900)(B5( I ) ,B6( I) ,  1=1,1000) 

CALL  MSEBET  (B5.B6, COUNT, XMSEPW) 

READ(51 ,900)(B7( I) ,B8( I ) ,  1=1,1000) 

CALL  MSEBET  (B7.B8, COUNT, XMSEBM) 

READ(61 ,900)(B9( I).BIO(I),  1=1,1000) 

CALL  MSEBET  (B9, BIO, COUNT, XMSEOLS) 

FORMAT  ( FI 1 .  8,2X,F11. 8) 

WRITE(6 ,*) 1 MSEDW1 
WRITE(6 ,*)XMSEDW 

WRITE(6,*) ' MSETN 1 
WRITE(6,*)XMSETN 

WRITE(6 ,*) 'MSEPW' 

WRITE(6,*)XMSEPW 

WRITE(6 ,*) 1 MSEBM1 
WRITE(6,*)XMSEBM 
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